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1 INTRODUCTION. 

The peculiar interest of magic squares and all lusus numerorum in general lies in the fact that 
they possess the charm of mystery. They appear to betray some hidden intelligence which by a 
preconceived plan produces the impression of intentional design, a phenomenon which finds its 
close analogue in nature. 
Paul Carus [3, p. vii] 

Magic squares have turned up throughout history, some in a mathematical context, others in 
philosophical or religious contexts. According to legend, the first magic square was discovered in 
China by an unknown mathematician sometime before the first century A.D. It was a magic square 
of order three thought to have appeared on the back of a turtle emerging from a river. Other magic 
squares surfaced at various places around the world in the centuries following their discovery. Some 
of the more interesting examples were recorded in Europe during the 1500s. Cornelius Agrippa 
wrote De Occulta Philosophia in 1510. In it he describes the spiritual powers of magic squares 
and produces some squares of orders from three up to nine. His work, although influential in the 
mathematical community, enjoyed only brief success, for the counter-reformation and the witch 
hunts of the Inquisition began soon thereafter: Agrippa himself was accused of being allied with 
the devil. Although this story seems outlandish now, we cannot ignore the strange mystical ties 
magic squares seem to have with the world and nature surrounding us, above and beyond their 
mathematical significance. 

Despite the fact that magic squares have been studied for a long time, they are still the subject 
of research projects. These include both mathematical-historical research, such as the discovery of 
unpublished magic squares of Benjamin Franklin [12], and pure mathematical research, much of 
which is connected with the algebraic and combinatorial geometry of polyhedra (see, for example, 
[1], [4], and [11]). Aside from mathematical research, magic squares naturally continue to be an 
excellent source of topics for "popular" mathematics books (see, for example, [3] or [13]). In this 
paper we explore counting functions that are associated with magic squares. 

We define a semi-magic square to be a square matrix whose entries are nonnegative integers and 
whose rows and columns (called lines in this setting) sum to the same number. A magic square is 
a semi-magic square whose main diagonals also add up to the line sum. A symmetric magic square 
is a magic square that is a symmetric matrix. A pandiagonal magic square is a semi-magic square 
whose diagonals parallel to the main diagonal from the upper left to the lower right, wrapped 
around (i.e., continued to a duplicate square placed to the left or right of the given one), add up 
to the line sum. Figure 1 illustrates our various definitions. We caution the reader about clashing 
definitions in the literature. For example, some people would reserve the term "magic square" for 
what we will call a traditional magic square, a magic square of order n whose entries are the integers 
1,2,. ..,n 2 . 

Our goal is to count these various types of squares. In the traditional case, this is in some 
sense not very interesting: for each order there is a fixed number of traditional magic squares. For 
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Figure 1: A semi-magic, a magic, and a symmetric pandiagonal magic square. 

example, there are 880 traditional 4x4 magic squares. The situation becomes more interesting if 
we drop the condition of traditionality and study the number of magic squares as a function of the 
line sum. We denote the total number of semi-magic, magic, symmetric magic, and pandiagonal 
magic squares of order n and line sum t by H n (t), M n (t), S n (t), and P n {t), respectively. 

We illustrate these notions for the case n = 2, which is not very complicated: here a semi-magic 
square is determined once we know one entry (denoted by * in Figure 2); a magic square has to 
have identical entries in each coefficient. Hence 
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Figure 2: Semi- magic and magic squares for n = 2. 



H 2 {t)=t + l 



M 2 (t) = S 2 (t) = P 2 (t) 



1 if t is even, 
if t is odd. 



These easy results already hint at something: the counting function H n is of a different character 
than the functions M n , S n , and P n . 

The oldest nontrivial results on this subject were first published in 1915: Macmahon [9] proved 
that 

dW V 4 /V 2 / 1° otherwise. 

The first structural result for general n was proved in 1973 by Ehrhart [7] and Stanley [15]. 

Theorem 1 (Ehrhart, Stanley) . The function H n (t) is a polynomial in t of degree (n — l) 2 
that satisfies the identities 



H n (-n -t) = (-l) n - L H n (t) , 



H n (-l) = H n {-2) 



H n (-n + l) =0 . 



The fact that H n is a polynomial of degree (n — l) 2 was conjectured earlier by Anand, Dumir, and 
Gupta [2]. An elementary proof of Theorem 1 can be found in [14]. Stanley also proved analogous 
results for the counting functions for symmetric semi-magic squares. In this paper, we establish 
analogues of these theorems for the other counting functions mentioned earlier. 
A quasi-polynomial Q of degree d is an expression of the form 



Q(t)=c d (t) t d + ■ ■ ■ + Cl {t) t + co(t) 
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where cq, c±, . . . , C4 are periodic functions of t and ^ 0. The least common multiple of the 
periods of the c,- is called the period of Q. With this definition we can state the first of our two 
main theorems. 

Theorem 2 The functions M n (t) , S n (t), and P n (t) are quasi-polynomials int of degrees n 2 —2n—l, 
n 2 /2 — n/2 — 2, and n 2 — 3n + 2, respectively, that satisfy the identities 

M n (-n -t) = (-I)™" 1 M n (t) , M n (-1) = M n {-2) = ■■■ = M n (-n + 1) = , 

S n (-n-t) = (-l) n(n - 1)/2 S n (t) , 5„(-l)=5„(-2) = -.- = 5„(-n + l) = 0, (1) 

P n (-n -t) = P n (t) , P n {-1) = P n {-2) = ■■■ = P n (-n + 1) = . 

We can extend the counting function H n for semi-magic squares in the following way. Define 
a semi-magic hypercube (also called a quasi-magic hypercube) to be a <i-dimensional n x • • • x n 
array of n d nonnegative integers that sum to the same number t parallel to any axis; that is, if we 
denote the entries of the array by rrij 1 ___j d (1 < jk < n), then we require XX=i m h-u = * f° r an 
k = 1, . . . ,d. Again we count all such cubes in terms of d, n, and t; we denote the corresponding 
enumerating function by H d (t). So H 2 = H n , whose properties are stated in Theorem 1. Except 
for the case d = n = 3, which seems to have appeared first in [5], we could not find any other 
references. We prove the following result for general d. 

Theorem 3 The function H d (t) is a quasi-polynomial in t of degree (n — l) d that satisfies the 
identities 

H d n {-n -t) = (-I)"" 1 H d n {t) , H d (-1) = H d (-2) =■■■ = H d (-n + 1) = . 

We prove Theorem 2 in Sections 2 and 3. To be able to compute the counting functions M n , 
S n , and P n for specific n, we need the periods of these quasi-polynomials. We describe in Section 4 
methods for finding these periods and hence for the actual computation of M n , S n , and P n . Finally, 
we prove Theorem 3 in Section 5. All of our methods are based on the idea that one can interpret 
the various counting functions as enumerating integer points ("lattice points") in certain polytopes. 



2 SOME GEOMETRIC COMBINATORICS. 

A convex polytope V in R d is the convex hull of finitely many points in M. d . Alternatively (and this 
correspondence is nontrivial [16]), one can define V as the bounded intersection of affine halfspaces. 
If there is a hyperplane i7={x€M (i : a-x = 6} such that V H H consists of a single point, then 
this point is a vertex of V . A polytope is rational if all of its vertices have rational coordinates. 

Suppose that V is a convex rational polytope in M. d . For a positive integer t, let L-p{t) denote the 
number of lattice points in the dilated polytope tV = {tx : x € V} (see Figure 3). Similarly, define 
Lp(t) to be the number of lattice points in the (relative) interior of tV. Ehrhart, who initiated the 
study of the lattice point count in dilated polytopes, proved the following [6]. 

Theorem 4 (Ehrhart) . If V is a convex rational polytope, then the functions L-p(t) and Lp(t) 
are quasi-polynomials in t whose degree is the dimension of V and whose period divides the least 
common multiple of the denominators of the vertices of V . 

In particular, if V has integer vertices, then L-p and LZ, are polynomials. Ehrhart conjectured and 
partially proved the following reciprocity law, which was proved by Macdonald [8] (for the case that 
V has integer vertices) and McMullen [10] (for the case that V has arbitrary rational vertices). 
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Figure 3: Dilation of a polytope. 



Theorem 5 (Ehrhart-Macdonald reciprocity law) . L-p(-t) = (— 1 



i dim V 



L*(t) 



Theorems 4 and 5 mark the beginning of our journey towards a proof of Theorem 2. The 
connection to our counting functions M n , S n , and P n is the following: all the various magic square 
conditions are linear inequalities (the entries are nonnegative) and linear equalities (the entries in 
each line sum up to t) in the n 2 variables forming the entries of the square. In other words, what 
we are computing is the number of nonnegative integer solutions to the linear system 

A x = b , 



where x is a column vector in W 1 whose components are the entries of our square, A denotes the 
matrix determining the magic-sum conditions, and b is the column vector each of whose entries is 
t. We will use the convention that the entries of x are ordered by rows of the original square, and 
from left to right in each row, that is, the entries xjk of the square are being arranged as 



X = (xn,X 12 , . ■ . ,X ln ,X 2 l,X 2 2, ■ ■ ■ ,X 2n , ■ ■ -,X n l,X n 2, ■ ■ ■ , X nn ) 
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The matrix A is constructed accordingly. For example, if we study magic 3x3 squares, 



A = 
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The first row of Ax = b represents the first row sum in the square: x\\ + x\i + x\% = t; the 
second and third row (X21 + ^22 + ^23 = t and X31 + X32 + £33 = t) are the sums of the second and 
third row in the square. The fourth row of Ax = b represents the first column sum in the square: 
xu + X21 + X31 = t; and the next two rows take care of the second and third column in our square. 
Finally the last two rows of Ax = b represent the two diagonal constraints X\\ + X22 + £33 = t and 

Xl3 + ^22 + Z31 = t. 

Furthermore, by writing b = t 1, where 1 signifies a column vector each of whose entries is 
1, we can see that our counting function enumerates nonnegative integer solutions to the linear 
system A x = t 1, that is, it counts the lattice points in the t-dilate of the polytope 

p = j x = (xi,. . .,x n 2) 6 M'" 2 : x k > for k = 1,2, . . . ,re 2 , A x = lj , 

provided that we choose the matrix A according to the magic-sum conditions. Note that V is an 
intersection of half-spaces and hyperplanes, and is therefore convex. No matter which counting 
function the matrix A corresponds to, the entries of A are Os and Is. We obtain the vertices of 
V by converting some of the inequalities x k > to equalities. It is easy to conclude from this 
that the vertices of V are rational. Hence by Theorem 4, M n , S n , and P n are quasi-polynomials 
whose degrees are the dimensions of the corresponding polytopes. Geometrically, the "magic-sum 
variable" t is the dilation factor of these polytopes. 

The Ehrhart-Macdonald reciprocity law (Theorem 5) connects the lattice-point count in V to 
that of the interior of V . In our case, this interior (again, using any matrix A suitable for one of 
the counting functions) is described by 

V int = jx = ( Xl ,...,x n 2) el" 2 : x k > for k = 1,2,... ,n 2 , A x = 1 j . 

The lattice-point count is the same as before, with the difference that we now allow only positive 
integers as solutions to the linear system. This motivates us to define by M*(t), S*(t), and P^(i) 
the counting functions for magic squares, symmetric magic squares, and pandiagonal magic squares 
as before, but with the restriction that the entries be positive integers. Theorem 5 then asserts, for 
example, that 

M n {-t) = (_i) dc g(*^) M*(t) . (2) 

On the other hand, by its very definition M*(t) = for t = 1,2, ... ,n — 1. Hence we obtain 
M n (—t) = for such t. Also, since each row of the matrix A defined by some magic-sum conditions 
has exactly n Is and all other entries 0, it is not hard to conclude that M*(t) = M n (t — n). 
Combining this with the reciprocity law (2), we obtain 

M n (t) = (-l)^(M n ) Mn{ _ t _ n) 
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There are analogous statements for S n and P n . Once we verify the degree formulas of Theorem 2, 
(1) follows from the observation that 

(-l)" 2 " 2 "- 1 = (-l) n_1 , 

(_l)f™ 2 -§™-2 = (_!)n(n-l)/2 ? 
^_]^™ 2 — 3n+2 _ ^_-Qra(n— 3) _ 

3 PROOF OF THE DEGREE FORMULAS. 

Let's start with magic squares: by Theorem 4, the degree of M n is n 2 (the number of places we 
have to fill) minus the number of linearly independent constraints. In other words, we have to find 
the rank of the n 2 x (2n + 2) matrix 
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Here we only show the entries that are 1; every other entry is 0. Similarly as in our example, the 
first n rows of A represent the n column constraints of our square (x\\ + • • • + x n \ = t, . . . , x\ n + 
■ • • + %nn = t), the next n rows represent the row contraints of the square, and the last two rows 
take care of the diagonal constraints. 

The sum of the first n rows of A is the same as the sum of the next n rows, so one of the first 
2ra rows is redundant and we can eliminate it; we choose the (n + l)th. Furthermore, we can add 
the difference of the first and (n + 2)th row to the (2n + l)th and subtract the nth row from the 
(2n + 2)th. These operations yield the n 2 x (2n + 1) matrix 
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These represent 2n + 1 linearly independent restrictions on our magic square. The polytope corre- 
sponding to M n therefore has dimension n 2 — 2n — 1: it lies in an affine subspace of that dimension, 
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and the point (1/rt, . . . , 1/n) is an interior point of the polytope. Hence, by Theorem 4, the degree 
of M n is n 2 - 2n - 1. 

The case of symmetric magic squares is simpler: here we have Y^k=i k = n ( n + places to 
fill, and the row-sum condition is equivalent to the column-sum condition. Hence there are n + 2 
constraints, which are easily identified as linearly independent. The dimension of the corresponding 
polytope, and the degree of S n , is therefore n 2 /2 — n/2 — 2. 

Finally, we discuss pandiagonal magic squares. Again we have n 2 places to fill; this time the 
constraints are represented by the n 2 x 3n matrix 
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Here the first n rows represent the column constraints in our square, the next n rows the pandiagonal 
constraints, and the last n rows the row constraints. 

As in the first case, the sum of the first n rows equals the sum of the next n rows, as well as 
the sum of the last n rows. We can therefore eliminate the (2n)th and (2n + l)th rows to get the 
n 2 x (3n — 2) matrix 

\ 
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In this new matrix, we can now replace the (n + l)th row by the difference of the (n + l)th and 
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first rows, the (n + 2)th row by the difference of the (n + 2)th and second rows, and so on: 
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Finally, we can replace the (n + 2)th row by the sum of the (n + l)th and (n + 2)th rows, the 
(n + 3)th row by the sum of the (n + 2)th and (n + 3)th rows, etc. Subtracting the (2n)th row 
from the (n + l)th row gives a matrix of full rank, that is, rank 3n — 2. The dimension of the 
corresponding polytope, which is the degree of P n , is thus n 2 — 3n + 2. 
This finishes the proof of Theorem 2. 

4 COMPUTATIONS. 

To interpolate a quasi-polynomial of degree d and period p, we need to compute p(d + 1) values. 
The periods of M n , S n , and P n are not as simple to derive as their degrees. What we can do, 
however, is to compute for fixed n the vertices of the respective polytopes, whose denominators 
give the periods of the quasi-polynomials (Theorem 4). This is easy for very small n but gets 
complicated very quickly. For example, we can practically find by hand that the vertices of the 
polytope corresponding to S3 are (2/3,0,1/3,1/3,2/3,0) and (0,2/3,1/3,1/3,0,2/3), whereas the 
polytope corresponding to £5 has seventy-four vertices. To make computational matters worse, the 
least common multiple of the denominators of these vertices is 60; one would have to do a lot of 
interpolation to obtain S5. 

The reciprocity laws in Theorem 2 essentially halve the number of computations for each quasi- 
polynomial. Nevertheless, the task of computing our quasi-polynomials becomes impractical with- 
out further tricks or large computing power. The following table contains data about the vertices 
of the polytopes corresponding to M n , S n , and P n for n = 3 and 4. It was produced using Maple. 



Polytope corresponding to 


Number of vertices 


l.c.m. of denominators 
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With this information, it is now easy (that is, easy with the aid of a computer) to interpolate 
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each quasi-polynomial. Here are the results: 
M 3 (t) = 



§i 2 + §t + l 







if 3\t, 
otherwise; 



s 3 (t) 
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if t = mod 4, 
if t = 2 mod 4, 
if + is odd; 



P 4 (t) 



ISO* 6 + ifo^ + l* 4 + * 3 + Is* 2 + T5* + 1 if * is even ' 



if £ is odd. 



5 SEMI-MAGIC HYPERCUBES. 

We now prove Theorem 3. The fact that H d is a quasi-polynomial, the reciprocity law for this 
function and H d *, and the location of special zeros for H d follow in exactly the same way as the 
respective statements in Theorem 2. As in that case, the remaining task is to find the degree of H d , 
that is, the dimension of the corresponding polytope. Again we have only to find the dimension of 
the affine subspace of W l in which this polytope lives. We do so by counting linearly independent 
constraint equations. 

To this end, let us write the coordinates of a point in the polytope corresponding to H d as 
c(ai, . . . , ad) with 1 < a,- < n. These are real numbers that satisfy the constraints 

n 

c(ai, . . . ,a d ) > , ^c(ai, . . . , o fc _i,j, a k+1 , . . . , a d ) = 1 (1 < k < d) . 

3=1 

Once the (re— l) d coordinates c(ai, . . . , a^) with 1 < a,j < re— 1 are chosen, the remaining coordinates 
are clearly determined by the foregoing conditions. Therefore, there are at least n d — (n — l) d linearly 
independent constraint equations. 

On the other hand, every constraint involves at least one of the n d — (n — l) d coordinates 
c(ai, . . . , ad) for which at least one dj equals n. Consider, say, the coordinate c(a\, . . . , a&, n,...,n) 
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in which < k < d and a±, . . . ,af- < n. We compute: 

n-l 

c(ai,...,ak,n,...,n) = 1 - ^ c(ai, . . . , a fc , n, . . . , n, j d ) 

n— 1 / n— 1 

= c ( a i'---' a fc' n '---' n 'i<i-i'id) 

id=i \ id-i=i 
= l-(n-l)+ ^ c(ai, . . . ,a fc ,n, . . . ,n,j d _i, j d ) 

l<jd-i ,3d<n-l 

= ... = 1 - (n - 1) + (n - l) 2 + (-l^'-V - 

+ (_1)<*-* c( ai ,...,a k ,j k+1 ,...,j d ) 

i<ife+iv,id<™-i 

1 - (1 -n) d ~ k . ,. d _ k „ 
= + (-1) >^ c(oi,...,o fc ,jfc+i,...,j ( i . 

We could have changed the order of summation, that is, used the d — k constraints that are in force 
here in a different order. However, we would always end up with the same expression. This implies 
that all constraints involving c(a\, . . . , a,}-, n, . . . , n) with < k < d and a±, . . . , a k < n, but not 
involving any coordinate with a,j = 1 for some j with j < k, are equivalent to the one constraint 



c(ai,...,a k ,n, ...,n) = h (-1) > c(oi, . . . , a fe , j fc +i, • • • , 3d) ■ 

n 



l<jfe+i,---Jd<n-l 



Therefore there are at most n d — (n — l) d linearly independent constraint equations. 
Accordingly, the dimension of the polytope and the degree of H d are (n — l) d . 

6 CLOSING REMARKS. 

One big open problem is to determine the periods of our counting functions. The evidence gained 
from our data seems to suggest that the periods increase in some fashion with n. We believe the 
following is true: 

Conjecture 1 . The functions M n , S n , and P n are not polynomials when n > 5. 

Results in [1] lend further credence to the validity of this conjecture. As for semi-magic hypercubes, 
Bona showed in [5] that iff really is a quasi-polynomial, not a polynomial. We challenge the reader 
to prove: 

Conjecture 2 . The function H d is not a polynomial when n, d > 3. 
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